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ABSTRACT 

In astrophysical regimes where the collisional excitation of hydrogen atoms is relevant, 
the cross sections for the interactions of hydrogen atoms with electrons and protons 
are necessary for calculating line profiles and intensities. In particular, at relative 
velocities exceeding ~ 1000 km s , collisional excitation by protons dominates over 
that by electrons. Surprisingly, the H-H + cross sections at these velocities do not exist 
for atomic levels of n 4, forcing researchers to utilize extrapolation via inaccurate 
scaling laws. In this study, we present a faster and improved algorithm for computing 
cross sections for the H-H + collisional system, including excitation and charge transfer 
to the n ^ 2 levels of the hydrogen atom. We develop a code named BDSCx which 
directly solves the Schrodinger equation with variable (but non-adaptive) resolution 
and utilizes a hybrid spatial-Fourier grid. Our novel hybrid grid reduces the number 
of grid points needed from ~ 4000n 6 (for a "brute force" , Cartesian grid) to ~ 2000n 4 
and speeds up the computation by a factor ~ 50 for calculations going up to n = 4. 
We present (I, m)-resolved results for charge-transfer and excitation final states for 
n = 2-4 and for projectile energies of 5-80 keV, as well as fitting functions for the 
cross sections. The ability to accurately compute H-H+ cross sections to n = 4 allows 
us to calculate the Balmer decrement, the ratio of Ha to H/3 line intensities. We find 
that the Balmer decrement starts to increase beyond its largely constant value of 2-3 
below 10 keV, reaching values of 4-5 at 5 keV, thus complicating its use as a diagnostic 
of dust extinction when fast (~ 1000 km s _1 ) shocks are impinging upon the ambient 
interstellar medium. 

Key words: atomic processes - ISM: lines and bands. 



1 INTRODUCTION 

The extreme densities (both low and high) and high tem- 
peratures inherent in many astrophysical processes allow for 
atomic interactions beyond the realm of terrestrial consider- 
ation or experiments. Hydrogen, the most abundant element 
in the Universe, signals its presence through the production 
of emission (or absorption) lines, either via the recombina- 
tion of protons with electrons or the collisional excitation of 
hydrogen atoms by electrons or protons. When the medium 
is tenuous (~ 1 cm -3 ) and the relative velocity of interac- 
tion is high (~ 1000 km s _1 ), recombination becomes slow 
and collisional excitation dominates. 

Fast, astrophysical shocks are examples where colli- 
sional excitation becomes important. In particular, a class 
of shocks known as "Balmer-dominated shocks" are driven 
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by astrophysical pistons (e.g., supernova remnants, pulsar 
wind nebulae, novae) impinging upon the ambient interstel- 
lar medium, producing hydrogen lines observed in the Ly- 
man and Balmer se ries accompanied by a dearth of metal 
lines fsee lHen J20ld for a review). Both the processes of exci- 
tation and charge transfer are inferred to be at work in these 
Balmer-dominated shocks, w hich may also be relev ant in 
young, high-redshift galaxies (|Heng fc Sunvaevl feoOS'). The 
line widths and intensities serve as diagnostics for the shock 
velocities and temperatures, which necessitates the knowl- 
edge of excitation cross sections to hydrogen levels n ^ 4, 
especially if one is interested in emission lines such as H/? 
as well as the Paschen and Brackett lines. The ratio of Ha 
to H/3 line intensities, known as the "Balmer decrement", 
further serves as a diagnostic for dust extinction du e to its 
insensitivity to electron temperature and density (|Drainel 
l201lf ) — it is of interest to establish if this insensitivity ex- 
tends to the regime occupied by fast, astrophysical shocks. 
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At relative velocities of ~ 1000 km s _1 or greater, 
the excitation of hydrogen atoms becomes dominated 
by the interactions w ith protons rather than electrons 
l|Heng fc McCravll2007h . A glaring and surprising gap in the 
existing literature, both from theory or experiment, is that 
cross sections for the reaction, 



H(ls) + H + -» H (n ^ 4) + H 4 



(1) 



are essentially non-existent at these velocities, correspond- 
ing to energies of about 5 keV. In the absence of these cross 
sections, some researchers have resorted to using approxi- 
mate scaling laws such as 



(2) 



where no < n, to extrapolate for cross sections with n ^ 4 
using available ones with no ^ 3. This scaling law can be 
derived using the Born approximation and is only approx- 
imately val id at high velocities v ^> ac, or E co i ^> 25 keV 
(Salin 1988). Here by velocity we mean relative velocity be- 
tween colliding particles, E co i is the kinetic energy of a hy- 
drogen atom moving towards the proton at rest. Besides the 
inaccuracy associated with extrapolation, it also leaves open 
the question of how to obtain cross sections for levels with I 
values which do not exist for no ^ 3 (e.g., 4f). 

Initial attempts to calculate and measure cross 
sections of hydrogen collis io ns date back t o the 
1960s (IStebbings et all Il965l; IWilets fc Gallaherl [iliol 



iRvding. Wittkower. fc Gilbodvl Il966l ; iBavfieldl Il969h . On 



the theoretical side, the success of these efforts has been 
hindered by the high computational cost of numerical 
simulations. On the experimental side, it has been limited 
by stringent requirements of creating high vacuum states 
and the high costs of preparing and characterizing atomic 
hydrogen targets. 

In calculating cross sections for high-ni proton- 
hydrogen collisions, it is important to consider several dis- 
tinct cases. At low velocities (v <C ac), collisions between 
a hydrogen atom and a proton lead to large deflections of 
the colliding particles and, as in the case of an unbound H% 
molecule, the electron wave function deforms adiabatically 
during the collisional time. In such situations, the initial 
configuration of the system is considerably modified during 
the collision and the process must be treated in a way which 
reflects the interplay between various quantum states of the 
electronic wave function. This problem has been addressed 
with the close-coupling approximation which assumes that, 
during the atomic collision, the electron wave function tran- 
sitions between a certain numbe r of configurations which 
form the "basis set" of functions (|Fritsch fc Linl fl991'). The 
dominant outcome of the low energy collisions is charge ex- 
change/transfer between the colliding particles. 

At high energies (v 2> ac), which includes the 
relativistic regime, colliding particles follow undeflected, 
straight-line trajectories. This case is well described by the 
Born approximation in which the incoming proton is seen 
as a small pertu rbation of the electronic wave function 
jBates fc G riffmg 1953J). The dominant outcome of the colli- 
sions in the relativistic regime is the excitation or ionization 
of the hydrogen atom. 

The third case corresponds to intermediate veloci- 
ties (v ~ ac), where collisional times are of the or- 
der of the atomic timescale and therefore a perturbative 



treatment of the problem becomes invalid. This is pre- 
cisely the regime relevant to Balmer-dominated shocks. 
In this regime, the behaviour of the electron wave func- 
tion is more complicated than at low energies or rela- 
tivistic energies. At intermediate energies (E ~ 10 keV), 
there is no dominant outcome for a collision: charge trans- 
fer, collisional excitation and ionization are all important 
and interconnected. There is no clear intuitive picture of 
the mechanism for populating various electron quantum 
states and therefore the use of the close-coupling approx- 
imation is challenging and requires development of mul- 
tiple basis sets and extensive convergence tests. Analysis 
of the proton-hydrogen collisions at intermediate energies 
in the cl ose-coupling approx i mation is an active area of 
research dFritsch 
iKuang fc Linlll996; 



Linlll99ll; iFord. Reading, fc Hall 
McLaughlin. Winter, fc McCann 



1993; 



1997 



Marti n 1999; Toshima 1999; Winteijl2009l ; ICrothers fc Dubd 
1 19921 ; iBrown fc Crothers! 1 19961 ). However, the accuracy 
of the obtained results is still not fully determined as 
convergence tests of these methods are extremely hard, 
especially if extended over a large range o f energies 
l|Ford. Reading, fc Hallll 19931 : IKuang fc Linlll996h . 

An alternative method to addressing the problem of 
proton-hydrogen collisions at intermediate energies is via 
the direct solution of t he Schrodinger differential equa- 
tion on a numer i cal gr i d (iMaruhn-Reswani, Grun. fc Scheid 
1979 1; iBottcherl Il982l ; iKulander. Sandhva Devi, fc Kooniii 
1982). In fact, the grid-based method can be thought of 



as a finite basis set method with one basis function for 
each point on the grid. Over t he past decade, several 
groups have taken this approach dKoIakowska et alj 1 19981 ; 
iKolakowska. Pindzola. fc Schultj Il999h . producing results 
for energies ranging from 10 to 100 keV. It is easier to test 
the convergence properties of these numerical grid methods, 
but the price to pay is that they are are notoriously com- 
putationally demanding due to the long-range nature of the 
Coulomb electrostatic force. This is especially true if one 
needs to accurately represent states of high n; for this rea- 
son, previous results were limited to n Sj 3. 

In most cases, the results obtained in previous studies 
measure and/or calculate cross sections of hydrogen colli- 
sions in the velocity ran ge ~ 100-1 000 km s" 1 only with a 
precision ~ 10-30% (see lHeng||2010l and references therein). 
There is also a substantial disagreement betwe en experimen- 
tal results and th eoretical calculations (e.g., IWinterl [20091 ; 
ISidkv fc Linll200ll ). In the case of final states with n > 3, 
robust theoretical or experimental cross-sections in the en- 
ergy range relevant to Balmer-dominated shocks studies do 
not exist at all. 

The objectives of the present study may be concisely 
stated as follows: 

• To introduce a novel hybrid grid (Figure [TJ for the di- 
rect solution of the Schrodinger equation; 

• To demonstrate that the use of this grid reduces the 
number of grid points needed from ~ 4000n 6 (for a "brute 
force" , Cartesian grid) to ~ 2000n 4 , which corresponds to a 
gain in the speed of computation by a factor ~ 50 for n ^ 4 
case; 

• To provide cross sections for excitation and charge 
transfer reactions, in the H-H + collisional system up to 
n = 4, at energies of 5-80 keV; 
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• To provide fitting functions for these cross sections so 
as to enable their (convenient) use by astrophysicists and 
astronomers; 

• To quantify the error associated with using the scaling 
law from equation ((2}; 

• To calculate the Balmer decrement in the regime where 
fast (~ 1000 km s _1 ) astrophysical shocks are impinging 
upon ambient interstellar medium (~ 1 cm -3 ). 

The rest of the paper is organized as follows. In Sec- 
tion [2] we provide a detailed description of the theoretical 
model behind our analysis and the major constraints driv- 
ing the development of the code for high-ni cross section 
calculations. Section [3] describes the code developed for our 
calculations and shows the results of extensive consistency 
tests. We discuss capabilities and limitations of the code and 
provide guidance on how this code can be expanded and used 
by other groups. In Section [3] as well as in Appendixes A 
and B, we provide results of our cross section calculations 
and compare our results with earlier studies. In Section [5] 
we briefly discuss astrophysical applications of the obtained 
cross sections with a specific focus on Balmer-dominated 
shocks. Our results are summarized in Section [6] 



2 COMPUTATIONAL MODEL FOR CROSS 
SECTION CALCULATIONS 

2.1 Initial setup for precise cross section 
calculations 

Our objective is to determine the cross sections for reactions 
of the form: 

H A (ls) + H+^X. (3) 

The hydrogen nuclei are assumed to be very massive so that 
their motion can be treated classically and one only has to 
follow the evolution of the electron wave function in the 
potential created by the two nuclei. Except for extremely 
small impact parameter b ^ rrip 1 v~ 2 , all potential energies 
are negligible compared to the nuclear kinetic energy, so 
we may use an undeflected (straight line, constant velocity) 
trajectory for the nuclei. In this limit, one can distinguish 
the two nuclei — hence their description as Ha and Hb- 

The initial electronic state is that of the Is orbital of 
atom A, i.e., [Isa). The final states X under consideration 
correspond to (i) no reaction (|1sa)); (h) excitation (\nlm,A), 
n ^ 2); (iii) charge transfer (\nlrriB)); and (iv) ionization 
(everything else). All of these are of interest, even for high n 
levels, for example for the H/3 lines the upper level is n — 4 
and for the Bra lines the upper level is n = 5. 

We choose a coordinate system such that the relative 
velocity points along the z-axis, V — Ve z , and the nuclear 
separation vector lies in the iz-plane. The relative separa- 
tion is 

va - r B = be x + Vte z . (4) 

Note that the electron wave function is always symmetric 
under reflection across the rrz-plane. We will choose the ori- 
gin of the coordinate system in the x-direction such that 
xa = and is = — \b. The choice of origin in the in- 
direction will be discussed later. 



The cross section to produce a particular final state X 
is given by 

r°° i „ |2 

<j x = lim / 2tt6 \{X\S(-T,T)\1sa)\ db, (5) 

where S{ti,tt) is the time evolution operator from time 
ti to tf. The S- matrix element can in principle be ob- 
tained by evolving |1sa) forward in time, \X) backward in 
time, or some combination of both. For example, one could 
evolve both states to t = 0, that is one could factor S as 
S(-T,T) = S(0,T)S(-T,Q) and have S(0, T) back-operate 
on (X\. The results in the present paper are based on evolv- 
ing | Isa) forward as this is the most efficient way to generate 
cross sections for large numbers of final states. 

Computation of the matrix elements requires us to solve 
the Schrodinger equation. In this work, we focus on grid 
methods because they allow for easier convergence tests and 
can be applied over a wide range of energies without signifi- 
cant modifications. Numerical grid methods also allow direct 
visualization of wave function evolution during the collision. 

Throughout this paper we use atomic units for all quan- 
tities, i.e., energy in hartrees, length in Bohr radii, velocity 
in units of ac w 2190 km s _1 , and mass in electron masses. 
Conversions between the Systeme Internationale (SI) units 
and atomic units are provided in Table [T] 

2.2 Grid choice 

The most obvious way to implement a grid method is to 
choose a spacing Aa; and a box size L. The number of grid 
points is then N ~ (L/Ax) 3 . Unfortunately this will be 
computationally prohibitive: if we want to consider a highly 
excited state of hydrogen nl, then the grid must go out to 
at least a radius of 2n 2 , and preferably much more, so L > 
4n 2 . On the other hand, to resolve the Is state properly a 
fine spacing (below ~ 0.2) is needed. This leads us to the 
conclusion that we need TV ~ 4000n 6 grid points, which is 
prohibitive for states above n = 4. 

Clearly, we will need a type of grid that puts resolu- 
tion where we need it: high resolution near the protons, and 
more modest resolution far away. To avoid the complexity 
of developing an adaptive code, we will insist on high res- 
olution near the trajectories of the proton and lower res- 
olution elsewhere. This immediately suggests developing a 
generalized cylindrical coordinate system, i.e., introducing a 
mapping (u, v) o (x, y) and using as our fundamental coor- 
dinates (u,v,z) instead of (x,y,z). A constant grid spacing 
Ait = Af can then correspond to a variable spacing in the 
xy-plane, in accordance with the Jacobian of the transfor- 
mation. Note that high resolution (several grid points per 
Bohr radius) is also required in the region in between the 
protons in order to correctly model the Isa Isb tunnel- 
ing that is primarily responsible for charge transfer at low 
and intermediate velocities. 

A minimal criterion for such a grid is that it should be 
able to adequately sample all bound wave functions with 
several grid points per cycle. The momentum of a bound 
wave function can be as large as \/2/p, where p is the min- 
imum separation from the nucleus. Therefore the grid spac- 
ing should be at most ~ \J p/2 (and preferably better). This 
requirement could be relaxed if p ^ 1, where the classical 
intuition concerning the "maximum momentum of a bound 
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Dimension 


Name 


Expression 


Value in SI 


Length 


Bohr radius 


ao = H/ (m e ca) 


5.29 X 1CT 11 m 


Energy 


Hartree 


Eh = a 2 m e c 2 


4.36 x 1CT 18 J 


Velocity 




ac 


2.19 X 10 6 m s" 1 


Electric field 




E h /(ea ) 


5.14 x 10 11 V m" 1 



Table 1. Connection between the atomic units and the SI units. 



electron" is invalid; in this regime instead all that is required 
is to have at least a few grid points per Bohr radius. 

After considering and rejecting several other choice^] 
we decided on the coordinate system 



1 + 



y = v\ 1 + 



(6) 



where u s is a parameter. The Jacobian is 



dx _ 1 + 2u 2 /u 2 s 
du ~ ^/l + u 2 /u 2 



1 

2u/u s 



\u\, \x\ <§C Us 
\u\, \x\ 3> U s 



and 



dy _ l + v 2 /2 
dv ~ ^/l + v 2 /4 



< i 
> i. 



(7) 



(8) 



This satisfies our resolution criteria if Au is no more than a 
few tenths, and u a is at least as large as ~ max(l, b). 

We illustrate the prosed approach in Figure[T]by using a 
sample case of A u = A z = 0.2 and comparing our proposed 
grid to a "brute force" Cartesian grid with = A z = 0.2. 
This figure shows every third point in both a;-direction and 
z-direction, and clearly illustrates benefits of putting high 
resolution in the region surrounding colliding particles while 
reducing the resolution far from the collision region. The 
resolution in the physical x-space between colliding particles 
is close to A x « A u — 0.2, while at larger separation the 
spacing between the grid points in z-direction significantly 
increases reaching A x ~ 1.6 at x ~ 16. 

By using a grid spacing of e.g. 0.2 in (u,v,z) space, 
and noting that we only have to go out to maximum values 
of u ma x ~ y/2u 3 n, Umax ~ v2n, the number of grid points 
necessary would be 



N ~ i(10\/2^rn)(10\/2n)(20n 2 ) 



2000uy 2 n 4 . 



(9) 



This is much more manageable than 4000n 6 as found earlier 
but is probably still too large for cases above n > 4. 

We can make yet another improvement by Fourier- 
transforming the z-direction. We suppose that we have a 
box of size L z in the z-direction and at each point in the 
(u, i>)-plane we do a Fourier transform, 



^(x, y, z) 



(10) 



where q is an integer. The maximum required value of q is 
now p z ,uiaxL z /2n. Here p z , m ax can be set to a large value 
(p z ,max ~ 20) within a few Bohr radii of the atoms (i.e., 
|u|,|u| less than a few), but a smaller value (p z , m ax ~ 4) 

1 For example, parabolic cylinder coordinates would have pro- 
vided the desired resolution for head-on collisions, but at large 
impact parameter would have difficulty providing resolution at 
the locations of both protons. 



at large separations. We refer to this as a q-drop procedure 
and it effectively reduces the number of grid points involved 
in the computation by removing Fourier modes that have 
very little contribution to the final result. For the results 
presented in this paper (cross sections with n = 4 final prin- 
cipal quantum number), it is computationally feasible to do 
the calculations at a single q m ax- However the code we have 
written, BDSCx, supports the use of two values of q mxc in dif- 
ferent regions of the (u, t>)-plane. The latter capability will 
be required for cross sections to the n ^ 5 levels of hydro- 
gen. We are currently tuning the parameters of the code 
for this purpose, but this paper presents the n sj 4 results 
with a single (/max in order to make these available to the 
community in a more expedient way. 

A peculiar property of this setup is that since we can 
only track momenta out to p z , m ax it actually matters which 
nucleus we take as moving and which as fixed. One would 
expect that the best results would be obtained by taking the 
H atom as fixed and the H + ion as moving, but for charge 
transfer reactions the definitions of "fixed" and "moving" 
change. We prefer to handle charge transfer by applying a 
boost operation by an amount V (i.e. increment the values 
of q by VL z /2tt) to the wave function of the electron. An 
alternative would be to evolve final states \nlms.) backward 
from t = T — > 0, apply a boost, and compute an inner 
product with the forward-evolved |1sa) state. For Hermitian 
discretized Hamiltonians these two methods are equivalent. 

To properly apply boost operation we consider two ref- 
erence frames: first, a stationary frame K with our initial 
coordinates (it, v, z) and the second is the restframe of a 
moving atom K', which moves with the velocity V in the 
+z direction. The transformation is 



*(r,t) = tf'(r- Vt,t)e' 



i(V-r-V 2 t/2) 



(11) 



2.3 Operators 

Here we restrict ourselves to orthogonal coordinate systems, 
i.e., where Vu ■ Vu = 0. 

We next need a method to compute inner products, 
and Hermitian discretizations of the kinetic and potential 
operators. The inner product is simply 



(12) 



where g uu = |Vtt| 2 , and g vv = \X7v\' 2 . 

The kinetic energy operator is the sum of the opera- 
tors along the x, y, and z axes. The z-operator is trivial, 
being simply a multiplication by 2ir 2 q 2 /L z . The x and y 
operators are trickier; fortunately, they commute with the 
Fourier transform in the z direction so we may implement 
them independently on each g-slice of the wave function. 
We recall that the x and y components of the kinetic energy 
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-5 5 

z-direction 

Figure 1. An example of a grid in xz-plane used in our anal- 
ysis (upper panel) with every third point plotted in both x- 
and z-direction, so that the aspect ratio is illustrated correctly. 
We also show positions of the moving hydrogen atom (thick dot 
with the vector) and stationary proton at Z = 0. Distances be- 
tween grid points in x-direction correspond to equidistant inter- 
vals (A u = 0.2 in u-space), so that in the region close to the 
interacting particles A x fa A u = 0.2, while at large separations 
spacing in a'-direction is significantly increased reaching, for ex- 
ample, Aj ~ 1.6 at x ~ 16. The lower panel shows a "brute force" 
Cartesian grid with A x = A z = 0.2 (again only every third point 
is plotted) and allows to clearly see the advantage of our grid 
choice. 



operator can be written as: 

_Ly^ f f uu dip* dip vv dip* dip 
2L Z J \ 9 du du 9 dv dv 



du dv 



(13) 



Thus, we see that the u and v parts of the kinetic energy 
operator are simply additive, i.e., we can write T — T u + 
T v + T z . The ii-part can be re-cast by discretizing the partial 
derivative as 



dip _ ip(u + Au/2) - tp(u - Au/2) 
du Au 



(14) 



Note that this partial derivative is measured not on the grid 
points, but halfway in between (i.e., Au/2 to the "right" of 
each grid point, or alternatively along each grid segment). 
Then we may write {ip\T u \ip) as a sum over such grid seg- 
ments. This gives an approximation to T u , 



t„v(«) = 



2Au 2 



g vv (v) 



tp(u + Au) — ip(u)] 



+ 



Au \ 
2 i 



9 vv (v) 



[tp{u) — ip(u — Alt)] 



(15) 



A simple calculation shows that with the discretized inner 
product of Eq. (|12p . this kinetic energy operator is exactly 
Hermitian. A similar equation holds for T v . Off-grid points 
are assumed to have ip = 0, corresponding to Dirichlet 
boundary conditions. We note, that versions of this opera- 
tor with higher order accuracy can be constructed by using 
more than 2 points in the derivative described in Eq. (|14p . 
An intelligent boundary condition would have to be chosen 
at the endpoints; however, only the ionized electrons will 
reach the boundary and they will reflect off. The higher or- 
der derivatives can be used to reduce spurious oscillations 
between ns o np nd states. We relegate the study of the 
use of the higher order derivatives to a future work. 

The potential energy operator is local in 3-dimensional 
position space, but not in (u, v, g)-space. The potential at 
any position is given by 



V(t) 



^(x-b/2) 2 +y2 + {z-z A ) 2 
1 

y/(x + b/2) 2 +y 2 + {z-l^ 



(16) 



The most efficient way to implement the potential operator 
is to FFT ip q (u, v) in the z-direction, multiply by V(x, y, z), 
and perform an inverse FFT. For N z points, this implies 
~ InNz operations per grid point, which is manageable es- 
pecially since in most cases N z is small because of the q-drop 
procedure. 

We note that this choice of potential term in the Hamil- 
tonian is associated with numerical difficulties due to its 
divergent nature near the proton. We eliminate these diffi- 
culties by capping the potential with a continuous function 
near the origin. We adopt 



V(r) 



4^(9-5^ 



r < Ro 
r > Ro, 



(17) 



which was chosen so that the volumetric integral vanishes: 
/ (14 a pp C d — Vtruc)d 3 r = 0. The advantage of this is that 
spurious features in the potential near the origin will re- 
sult in a spurious interaction Hamiltonian between any two 
states tpi and ip 2 given by / V , i^2(V cap p C d - K ruc )d 3 r, and 
hence setting the volumetric integral to zero should yield 
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improved behavior over, e.g., imposing a simple floor on V. 
The capping radius is chosen to be _Ro = 0.2 for the results 
shown here. The capping procedure was tested extensively 
and shown not to introduce any spurious deviations from 
the results with uncapped potential. Details of the testing 
procedures are discussed in Subsection 13.31 and illustrated 
in Fig. [4] For all relevant quantities calculated using our 
code the difference between capped and uncapped potential 
is << 1%. 



3 THE CODE 

3.1 Grid parameters 

To implement the proposed algorithm we developed a 
grid-based code BDSCx, which computes cross sections of 
hydrogen-proton collisions with hydrogen starting in | Is) 
state: 



| Is) = *ioo = -j=e 

V7T 



(18) 



We start and end the collision when the the two parti- 
cles A and B are separated by a distance sufficient to fully 
resolve wave functions of interest. Distance from the edges of 
the box is also determined by the requirement to properly 
resolve nlm states of interest at the beginning and at the 
end of the collision. Our choice of the box size parameters is 
guided by the charge distribution in the states of interest. In 
Figure O we plot the charge density r 2 R^ t (r) for hydrogen 
states with n — 1,2,3. Here, R n i(r) are normalized radial 
eigenfunctions of the hydrogen atom: 



Rni(r) = 



(r 



I 



4n 4 [(n + ;)!] 3 



jvj e _r / r ' 



n 



21+1 
n-l-l 



2r 



(19) 

where i^ i 1 _ 1 (2r/n) are the generalized Laguerre polynomi- 
als. 

We choose the length of the box in the z-direction by 
requiring charge density decrease of more than 2 orders of 
magnitude relative to its maximum value for the state of 
interest. For example, to properly resolve the Is state a par- 
ticle should be separated from the edge of the box by more 
than 5 Bohr radii, whereas for the n = 2we need more than 
15 Bohr radii of separation. These conservative resolution 
requirements lead to the size of the box in z direction given 
by L z = 5 + 15 x 4 = 65 in the case when only n — 2 states 
are of interest. 

We also use the same conservative requirements for the 
box size in the x (u) and y (v) directions. The size in the y 
direction is fixed as soon as we decide on the upper n state 
of interest, while the size in x direction also depends on the 
impact parameter 6 so that L x — L y + b. We further require 
high resolution near the particles so that Au ~ Av ~ Az ~ 
0.18. This resolution requirement was tested and found to 
converge with the difference between A ~ 0.18 and A ~ 
0.16 being less than 0.1 per cent. Several examples of input 
parameters required for accurate cross section results are 
provided in Tabled 

One of the important advantages of the grid approach 
chosen in our calculations is the ability to visualize the evolu- 
tion of the wave functions during the collision. In Figure 



we show the time evolution of electron's probability den- 
sity as the hydrogen atom moves past the stationary proton 
with impact parameter b — 2 and velovity V = 1 a.u. We 
clearly see that after the collision part of electron's proba- 
bility density is spread between the two atoms, indicating 
possible charge transfer during the impact. 



3.2 Collisions with large impact parameter 

At large values of the impact parameter b the number of 
required grid points in the u and v directions increases dra- 
matically. Fortunately, at large separations (b ^ 5) the cross 
section results are well described by the Born approximation 
if the collision energies are reasonably high (E co \ ^ 1 keV). 
To reduce the amount of numerical computations without 
sacrificing precision of the cross-section results we use the 
Born approximation to obtain transition probabilities above 
b = 5. It is easy to show that in this approximation the tran- 
sition probability into the X final state can be written as: 



\(X\S(-T,T)\ls}\ 2 = 



i(E x -E )t 



(X\W\ls) dt 



(20) 

where W is the perturbation Hamiltonian operator, which 
in our case is given by 



w = -- 



1 



R |R-r| 



(21) 



and R is the vector separating the two protons. If one ap- 
proximates the Hamiltonian as a dipole, W ~ r • R/|i?| 3 , 
as appropriate at large impact parameters, then an analytic 
solution for the transition probability is obtained. In the first 
order, probability of excitation into \X > can be written as: 



KX|S(-oo,oo)|ls)| 2 w 
j H {X\x\ls) + 



i I e 

' — CO 



-i(E t -E a )t/v 



(22) 



(b 2 + £ 2 ) 3 / 2 



For example, the first order probability of transition 
into 2p state becomes 



|<2p |S(-oo,oo)|ls)| s 



/ 32>/2Ar (fS ; 

I 81V 2 



(23) 



where Ko(r) is a Bessel function of the second kind. 

Note that the Born approximation only allows transi- 
tions into the \npA) final states (and to ionized final states). 

At collision energies of E co \ > 1 keV and impact pa- 
rameters of b > 5, the results produced by using the Born 
approximation differ from the results obtained by our code 
by less than 5 per cent, and because of the small value of 
the transition probabilities at large 6 values the difference 
introduced into cross section results is <C 1 per cent. 



3.3 Convergence and consistency tests 

To test the proper functionality of our code and fully under- 
stand the dependence of the results on the input parameters 
such as grid resolution, the size of the box within which we 
keep a large number of point in q-space, box size and others 
we ran a large number of tests most important of which are: 
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n 


b 


U s 




Ly(L V ) 


L z 


N u 


N v 


N z 


2 


1 


1 


31(8) 


30(11) 


65 


46 


62 


362 


2 


5 


5 


35(18) 


30(11) 


65 


100 


62 


362 


4 


1 


1 


101(15) 


100(20) 


205 


84 


112 


1140 


4 


5 


5 


105(32) 


100(20) 


205 


178 


112 


1140 


5 


1 


1 


143(18) 


142(24) 


285 


100 


134 


1584 


5 


5 


5 


147(38) 


142(24) 


285 


212 


134 


1584 



Table 2. Examples of the simulation parameters for collisions involving n = 2 and n = 4 states. 



(i) Time-dependent evolution of an electron in the field 
of a single proton at rest. In this case we would expect that 
the state of a system remains unchanged, i.e., if a system 
starts in | Is) state it remains there and the probabilities of 
transitions into states with n = 2 vanish within the numer- 
ical precision of a discretized Hamiltonian system. Figure [4] 
shows the result of system's evolution over a few dynam- 
ical times (Tdyn = 27rn 3 ). We clearly see that the system 
exhibits small oscillations associated with discretization of 
the Hamiltonian; however, the results remain close to the 
expected values for |(ls|S(T)|ls)| 2 and for the energy ex- 
pectation {Ho}. For example, |(ls|5'(T)|ls)| 2 w 1 over the 
full time period with the precision better than 0.05 per cent; 
energy expectation deviates from the theoretical prediction 
(H ) = -|n~ 2 by less than 0.5%. The probability of Is) 
to 1 2s) transition remains small and oscillates with the pe- 
riod T osc ~ 37r/4 determined by the energy splitting between 
n — 1 and n — 2 levels. Furthermore, the precision of the 
results increases as we increase the resolution of the grid, 
allowing us to achieve desired accuracy of the final results 
through convergence. 

(ii) The next test incorporates the motion of the hydrogen 
atom. In this case, the initial wave function of the electron 
is phase-shifted by the factor e lVz as discussed at the end of 
Sec. 12.21 The results of these tests are also shown in Figure^ 
We clearly see that the normalization of the jls) is properly 
conserved. We note, however, that the motion of a particle 
introduces small additional oscillations into the shapes of 
the curves. 

(iii) We tested convergence of our cross section results 
by running our code with three different resolution values: 
A = 0.175,0.18 and 0.22 at collision energy E co i = 80 keV. 
The results converge, with cross section difference between 



1.000 r 

0.500 - 

£ 0.100- 
g 0.050 - 

1s 




0.001 

10 20 30 40 50 

Time 



Figure 4. Deviation of |(ls|S(T)|ls)| from unity as a function of 
time in a test runs without the second particle. Thick curves show 
the evolution of a stationary hydrogen atom located at x = 0.5, 
y = and z = —20 for the cases with unmodified 1/r poten- 
tial (solid) and modified potential of Eq. 1171 (dashed). Regular 
curves show the evolution of a hydrogen moving in the positive 
z-direction with V = 1 with unmodified potential (solid) and 
capped potential (dashed). Dash-dotted curve shows deviation 
from the exact energy expectation {Ho) = — l/(2n 2 ) for station- 
ary hydrogen. 

the A = 0.175 and A = 0.18 cases being less than 3%. We 
also checked that increasing the box size at a fixed resolu- 
tion does not modify our results, which indicates that the 
results of the boundary interactions are minimal for the set 
of parameters used in our runs. 

It is important to note that, while our cross section re- 
sults converge within the limits of simulation accuracy, they 
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Figure 3. Contour plots showing time evolution of the electron probability density as hydrogen atom moves with V = 1 and 6 = 2 past 
the proton at rest. Time in the panels is given in atomic units. 



still contain inherent uncertainty associated with our dis- 
cretization procedure, finite resolution of the numerical grid, 
finite size of the collision box, approximation of straight line 
trajectory for colliding protons, reconstruction of the cross 
sections from discrete probability values on a grid of im- 
pact parameters b, and the use of the Born approximation 
at high values of impact parameters. Numerical errors as- 
sociated with these sources are extremely hard to quantify 
precisely, although for each individual source the error has 
been tested and minimized with the convergence procedure 
to be below 1%. 

In testing we also compared our predictions for 
low n (n ^ 3) cross sections with the results of pre- 
vious studies and showed that our results are fully 



consis tent with the results presented in [Kuang fc Lin! 
i 19961). iKolakowska. Pindzola. fe SchultzT <|l999h . 



and IWinterl (|2009l ). We note, however, that the major 
goal of this paper is to introduce the computational 
algorithm and the code that allow accurate cross section 
calculations at high values of n as well as illustrate the 
importance of these calculations in studying Balmer- 
dominated shocks. Therefore, detailed comparison of our 
results to the results of other studies is not performed here 
and will be presented in future papers. 



4 CROSS SECTION RESULTS 

In this section, we present results for the cross sections for 
excitation and charge transfer final states, and (where pos- 
sible) compare our results to previous computations. 

We obtain cross sections by integrating Eq. ([5]), with 
\X) — \nlm). The final state wave functions in the rest frame 
of the final atom (A or B) are given by 



\nlm) = *nim(r, 



R n i(r)Yr(8, 



(24) 



Here, R n i (r) are normalized radial eigenfunctions (Eq. I19|) 
and Yi m (8,4>) are spherical harmonics. If the final atom is 
moving, then a boost must be applied to the relevant wave 
function. 



The first results obtained by running our code corre- 
spond to collisions involving only n = 1 and n = 2 states. 
This problem is the least computationally demanding and 
it has a large amount of data produced by previous studies, 
allowing for further tests of the code. In Figure [5] we plot 
probabilities of excitation to 2s state and charge transfer 
into Is state in collisions with impact energy E co \ = 40 
ke V. We compare our re sults with the results obtained 
in IKolakowska et al.l ljl99ct ) and show good agreement over 
the full range of b with slight deviation at small values of 
impact parameter caused by higher resolution used in our 
runs. Our deviation in the case of charge transfer into Is 
state causes higher value of the overall cross section, which 
is consistent with the resu lts obtained by other groups (see 
for example I Winter! (|2009h and references therein). 

A significant advantage of using direct solution of the 
Schrodinger equation on a grid is the ability to see the con- 
vergence of the probability results. In Figure [6] we show 
an example of probability evolution for excitations into 2s 
states. This method clearly allows to see the convergence of 
the numerical calculation and allows direct comparison of 
runs with various box parameters. 

We further note that our cross section results are fully 
consistent with the results obtained by other groups for low 
n cases. For example, in Figure [7] we plot cross section for 
charge transfer into Is state for various collision energies 
and compare our re s ults with other theoretical calculation s 
(|Kuang fc Linll 19961 ; IKolakowska et al. I ll998l ; IWinteri r2009). 

In this study we have limited our cross section calcula- 
tion to n = 4, which allows us to calculate Balmer decrement 
discussed in the next section. Our results, showing the cross 
sections of charge transfer and excitations are provided in 
Tables [3] and [4] In Appendixes A and B, we also provide 
Chebyshev polynomial fits to our results for excitations and 
charge transfer into 3s, 3p, 3d, 4s, 4p, Ad, and 4/ states 
and plots of obtained cross se ction s compa r ed to the results 
from IKolakowska et al.l (119981 ) and I Winter! (|2009l ). We note, 
that the BDSCx code produces the cross-sections for transi- 
tions into states with various n, I, and m, allowing for stud- 
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Figure 5. Plot on the left shows probability of charge transfer into Is state as a function of impact paramet er for E co i = 40 keV, an d 
the plot on the right shows probability of excitation into 2s state. Diamonds show corresponding results from Kolakowska ct al 
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Figure 6. Probability of excitations into 2s states as a function 
of time for 2 different values of impact parameter b = 0.2 a.u. 
(solid curve) and b = 0.8 a.u. (dashed curve). 



ies of polarization-dependent signals. Calculation of higher 
n states, as well as more detailed analysis of the g-drop pro- 
cedure are relegated to the future work. 



5 ASTROPHYSICAL APPLICATIONS 

5.1 Errors Associated with Extrapolating from 
Cross Sections with Lower n-values 

In the absence of available cross sections, one is often forced 
to extrapolate from known cro ss sections with lower princi- 
pal quantum numbers n (e.g.. lHeng fc McCravl [200 7:1, thus 
generating errors in the calculation of line profiles and in- 
tensities which are unquantifiable. With the benefit of now 
being able to calculate excitation cross sections for the H- 
H + collisional system, we quantify the error associated with 
using the scaling law in equation ([2]) for obtaining cross sec- 
tions with n 3. In Figure [SJ we see that the errors associ- 
ated with extrapolating for n — 4 cross sections from no = 3 
are typically a factor ~ 2. The less necessary extrapolation 
of obtaining n = 3 cross sections from no = 2 (since data 
for n = 3 is available) results in errors of a factor ~ 5. Since 
these errors are non-negligible, they probably dominate any 



o 

b" 




100 



Figure 7. Excitation cross sections (in units of 10 — 18 cm 2 ) for 
charge transfer into Is state. Stars show results of our calcula- 
tions, diamonds correspond to Kolakowska (1998), squares show 
results from Winter (2009) and plus signs correspond to the re- 
sults from Kuang & Lin (1996). 



uncertainty associated with a numerical integration tech- 
nique used to compute line profiles and intensities. A direct 
calculation, such as the one we have performed in this study, 
is necessary in order to obtain accurate cross sections and 
in turn perform a spectral analysis of lines such as H/3 at 
energies > 5 keV. 



5.2 The Balmer Decrement 

The Balmer decrement generally refers to the ratio of dif- 
ferent lines in the Balmer series: Ha/H/3, H/3/H7, etc. It 
is somewhat insensitive to the electron temperature and 
(low) densities. For ex ample, Ha/H/? has a value of about 
2-3 (|Osterbrocklll96Sl ). unless collisional excitation of hy- 
drogen atoms b y electrons dominate in wh ich case its value 
is as high as 8 (| Adams fc Petrosianlll974T ). which only oc- 
curs at electron densities ~ 10 4 cm -3 or higher. The rel- 
ative insensitivity of the Balmer decrement to the atomic 
processes makes it an appropriate diagnostic for the pres- 
ence of dust extinction, since the bluer lines in the Balmer 
series (e.g., H/3) are subjected to increased extinction and 
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1.0 


0.69 


0.54 


0.35 


0.05 


0.08 


0.05 


0.007 


6 x 10~ 4 


30 fceV 


1.2 


1.15 


0.88 


0.59 


0.36 


0.06 


0.08 


0.05 


0.007 


8 x 10~ 4 


40 keV 


1.35 


1.5 


1.1 


0.6 


0.4 


0.07 


0.06 


0.04 


0.007 


8 x 10 -4 


60 fceV 


1.15 


1.6 


1.4 


0.5 


0.3 


0.09 


0.03 


0.02 


0.005 


0.001 


80 fceV 


1.1 


1.6 


1.5 


0.3 


0.2 


0.08 


0.015 


0.008 


0.003 


0.001 



Table 3. Cross section results for excitation transitions into n = 2, 3 levels. 



cr/V 


Is 


2s 


2po 


2p±i 


3s 


3po 


3p±i 


3d 


3d±i 


3d±2 


5 keV 


1092 


5.8 


2.1 


11.5 


0.25 


0.35 


0.60 


0.20 


1.0 


0.008 


7.5 fceV 


928 


12 


3.0 


12.9 


0.55 


0.67 


0.89 


0.25 


1.2 


0.03 


10 keV 


795 


18 


4.9 


13 


1.5 


0.90 


1.5 


0.30 


1.1 


0.05 


12.5 fceV 


695 


27.1 


6.6 


11.6 


3.2 


1.6 


1.6 


0.40 


0.80 


0.04 


15 keV 


593 


32.8 


7.8 


9.8 


4.8 


2.1 


1.6 


0.38 


0.52 


0.03 


20 fceV 


425 


39 


7.9 


6.5 


8.6 


2.5 


1.7 


0.35 


0.28 


0.03 


25 fceV 


309 


38.6 


6.8 


4.4 


8.9 


2.4 


1.0 


0.28 


0.12 


0.015 


30 keV 


224 


34.9 


5.7 


3.0 


8.8 


2.0 


0.80 


0.21 


0.07 


0.01 


40 fceV 


120 


22 


3.6 


1.5 


6.5 


1.3 


0.40 


0.10 


0.03 


0.005 


60 keV 


42 


8.6 


1.3 


0.40 


2.7 


0.50 


0.15 


0.03 


0.01 


0.002 


80 fceV 


17 


3.5 


0.49 


0.15 


1.1 


0.20 


0.05 


0.01 


0.004 


5 x lO" 4 




4s 


4p 


4p±i 


4d 


4d±i 


4d± 2 


4/o 


4/±i 


4/±2 


4/±3 


5 fceV 


0.12 


0.20 


0.10 


0.20 


0.30 


0.003 


0.05 


0.20 


0.001 


2 x 10" 5 


7.5 keV 


0.086 


0.29 


0.14 


0.07 


0.40 


0.008 


0.04 


0.17 


0.005 


8 x 10~ 5 


10 fceV 


0.64 


0.39 


0.34 


0.11 


0.50 


0.02 


0.04 


0.11 


0.009 


1 x 10" 4 


12.5 keV 


0.86 


0.70 


0.45 


0.22 


0.40 


0.016 


0.039 


0.08 


0.005 


1 x 10 -4 


15 fceV 
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0.94 
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0.29 


0.30 
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0.05 


0.003 


1 x 10" 4 


20 fceV 
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1.2 
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0.17 
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0.03 


0.02 


0.003 


6 x 10~ 4 


25 keV 


3.5 


1.2 


0.38 


0.24 


0.08 


0.008 


0.02 


0.005 


6 x 10~ 4 


5 x 10" 5 


30 fceV 


3.6 


1.0 


0.30 


0.20 


0.05 


0.005 


0.01 


0.002 


3 x 10~ 4 


6 x 10" 5 


40 keV 


2.7 


0.65 


0.17 


0.18 


0.02 


0.004 


0.005 


0.0015 


1 x 10~ 4 


1 x 10" 5 


60 fceV 


1.2 


0.23 


0.06 


0.02 


0.006 


9 x 10~ 4 


9 x 10" 4 


2 x 10~ 4 


3 x 10" & 


3 x 10" b 


80 fceV 


0.50 


0.10 


0.02 


0.007 


0.002 


0.0003 


2 x 10~ 4 


6 x 10~ b 


1 x 10~ b 


1 x 10~ b 



Table 4. Cross section results (in units of 10 18 cm 2 ) for charge transfer transitions into n = 2, 3 and 4 levels. 



thus the Balmer decrement attains a value larger than, e.g., 
2-3 in the case of Ha/H/3. A high value of the Balmer 
decrement may also be attained if the population of ex- 
cited hydrogen (e.g., 2s) is sufficient to cause self-ab sorption 
in the Balmer lines l|Capriotti|[T96i ; lNetzerlll975l) . As ex- 
amples, the Balmer decrement has been used as a di- 
agno stic in the study of supern ovae (e.g.. lAldering et al.l 
|2006| ). active galactic nuclei (e.g.. iDong et all 20081 ) and the 



Sloan Digital Sky Survey (SPSS) sam ple of galaxies (e.g., 
iGroves. Brinchmann fc Walchej[20Tll ). 

Since our work was originally motivated by the study of 
fast astrophysical shocks, we use them as an example in the 
calculations presented in this sub-section. In particular, we 
use the example of Balmer-dominated shocks, which are fast 
(~ 1000 kms" 1 ) shocks impinging upon tenuous media (~ 1 
cm" 3 ). As Balmer-dominated shocks are mostly observed 
around young (< 1000 years), Galactic supernova remnants, 
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Figure 8. Errors associated with the calculated versus extrapo- 
lated (via equation [2J cross sections as a function of the impact 
energy. Upper panel shows percentage errors for 3s (solid curve), 
3po (dashed curve), 3pi (dotted curve), and 4s (dot-dashed curve) 
excitations. Lower panel shows errors for 4po (solid curve), 4pi 
(dashed curve), 4do (dotted curve), 4di (dot-dashed curve), and 
4^2 (thick dot-dashed curve) transitions. 
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the production of photons via radiative recombination is 
unimportant since the recombination time is ~ 10 4 years. 
The calculation of line intensities then requires knowledge 
of how the various atomic levels of hydrogen are populated 
via collisions as well as their subsequent rates of radiative 
decay. For a strong shock (Mach number greatly exceeding 
unity), the relative velocity between the electrons/protons 
and hydrogen atoms is 8v = 3t> s /4, where v s denotes the 
shock velocity. Thus, the shock velocity can be related to 
the interaction energy E co i via 



1300 km s~ 



E, 



:ol \ 



1/2 



5 keV J 



(25) 



where ran is the mass of the hydrogen atom. The cross sec- 
tions presented in the present study are thus relevant to 
shocks with v s « 1300-5200 km s" 1 . 

Formally, evaluating the rate at which a given nl level 
is populate d by collisions requires the calculation of the rate 
coefficient jHeng fe McCravll2007T ), 



n 



Fi(vi) F 2 (v 2 ) Ava(Av) d 3 «i d 3 « 2 . (26) 



The preceding, six-dimensional integral is evaluated over all 
of the relative velocities (Aw) between the velocity distribu- 
tions of the atoms (Fi) and electrons/protons (F2), weighted 
by the relevant cross section (<r) of the atomic process being 
considered. 

To gain an intuition for the functional dependence of 
the Balmer decrement on the interaction energy E co i , it is 
sufficient to consider either individual pairs of particles or 
particles in Delta-function distributions. In this case, the 
rate coefficient reduces to 1Z = aAv. The Balmer decrement 
is then the sum of the cross sections for the collisional popu- 
lation of each nl level weighted by the appropriate branching 
ratio, 



Ha/H/3 



^(3s) + B 3 p,2s^(3p) + ^(3d) 



B 4s,2p CT ( 4s ) + B 4p,2s CT (4p) + B 4d)2p <r (4d) - 

(27) 

The branching ratio is simply the Einstein A-coefficient for a 
given transition normalized by the Einstein A-coefficients of 
all of the transitions allowed by the electric dipole selection 
rule. For example, B 3 p 2s = ^3p,2s/(^3p,2s + A 3p,ls)- 

Two extreme lim iting cases are typica ll y considered : 
Case A and Cas e B l|Baker fe Menzell Il938l ; ISeatonl li"960l ; 



Figure 9. The Balmer decrement (Ho/H/3) as a function of the 
impact energy for collisional excitation of hydrogen atoms by pro- 
tons. Solid line corresponds to Case B, while dashed line shows 


B 3p,2s R 
B 4s,2p R 


i 0.1183, 
t 0.5841, 


Case A. 


B 4p,2s R 


i 0.1191, 




B 4d,2p f 


« 0.7456 



I II "T 

Osterbrock 1989). Case B occurs when the neutral hydrogen 
column density is large enough for Lyman lines to be opti- 
cally thick. They undergo multiple scatterings and are even- 
tually degraded into a Balmer line and Lya or two-photon 
emission. Case A occurs when the neutral hydrogen column 
density is small enough for Lyman lines to be optically thin 
and s o freely escape the clou d. In a Case A scenario, we have 
(e.g., iHeng fc Sunvaevll2008l ), 



(28) 



In a Case B scenario, these branching ratios are essentially 



unity, since 3 S ~ 10 and B^p 3 S 



10" 
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In Figure [9] we show calculations of the Balmer decre- 
ment, for both Case A and B, using equation ([27} • For sim- 
plicity, we consider only the collisional excitation (and not 
charge transfer) of hydrogen atoms by protons; collisional 
excitation by electrons is sub-dominant at these energies. 
It is apparent that Ha/H/3 remains somewhat constant at 
values of 2-3 at high energies and starts to grow rapidly at 
energies below ~ 10 keV. At 5 keV, it reaches values of 4-5. 
If we use the values of 2-3 as a baseline, then this corre- 
sponds to the true dust extinction being over-estimated by 
AAy =1-3 (based on Ry = 3 .1 mo del and extinction curves 
from IWeingartner fc Draind l|200ll )L Thus, we caution the 
use of the Balmer decrement as a diagnostic for dust extinc- 
tion as it possesses some sensitivity to the atomic physics at 
energies < 10 keV. 



6 SUMMARY 

In this work, we introduced a new formalism for computing 
precise cross-sections for high-nZ proton-hydrogen collisions 
and developed a numerical code which implements our for- 
malism. We further used our code to obtain accurate cross 
sections for collisions between protons and hydrogen atoms 
which start in the ground state. Our computed cross sec- 
tions focused on the energy range of direct interest for the 
studies of Balmer-dominated shocks; as the observed spectra 
of these shocks improve in quality and precision, our cross 
sections are required for doing a detailed interpretation, e.g., 
to estimate the degree to which electron and proton temper- 
ature are equilibrated. 

The code, BDSCx, introduced in this paper has a large 
number of potential applications in atomic physics. Using 
the g-drop procedure and curvilinear coordinates proposed 
in our work enables relatively inexpensive calculations of 
charge transfer and excitation cross sections for proton- 
hydrogen collision with n < 7. In this paper we focused 
on the formalism and numerical implementation of the pro- 
posed cross section calculations, while detailed tests of the 
(/-drop procedure and the cross section results for n > 4 will 
be reported in future papers. 
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APPENDIX A: FITTING FUNCTIONS FOR 
COMPUTED CROSS SECTIONS 

In order to facilitate the broader use of our cross sections, we 
are providing fitting functions which allow a more straight- 
forward utilization of the obtained results. We fit the ob- 
tained cross sections for excitations and charge transfer into 
3s, 3p, 3d, 4s, 4p, 4d, and 4/ states with a series of Cheby- 
shev orthogonal polynomials. The fitting function is of the 
form: 

F{x;A)=exp(^+^ AiCi(x)\ (Al) 



where the coefficients A — Ai for ^ i ^ 7 are the fitting 
parameters. The quantities d are the Chebyshev orthogonal 
polynomials: 



Ci (x) = X, 


(A2) 


C 2 (x) = 2x 2 - 1, 


(A3) 


C 3 (x) = 4a; 3 - 3a;, 


(A4) 


d (x) = 8 (x A - x 2 ) + 1, 


(A5) 


C 5 (x) = 16x 5 - 20x 3 + 5x, 


(A6) 


C (i (x) = 32a; 6 - 48a; 4 + 18x 2 - 1, 


(A7) 


C 7 (x) = 64x 7 - 112x 5 + 56a; 3 - 7x. 


(A8) 


The fitting variable x is defined as 





^ _ In (Ecol/Emin) — In (Em^/Ecl) 
In (-Emax/-Bmin) 

where E co i is the relative energy between the proton and 
hydrogen atom; £J m j n = 5 keV and i? max = 80 keV are the 
respective minimum and maximum energies in our simula- 
tion. We use the Newton fitting algorithm which provides 
a local-optimal fit to the array of available cross section re- 
sults. 

The fitting parameters are provided in Tables IA1I 
and IA2I We note, that for purposes of precise analysis the 
use of actual cross section data points provided in Tables [3] 
and [3] and simple spline extrapolation will produce more ac- 
curate results. 



APPENDIX B: CROSS SECTION PLOTS 

This Appendix provides plots of cross sections for excitation 
and charge transfer transitions into 3s, 3p, 3d, 4s, 4p, Ad, 
and 4/ states. Our results are co mpared with the results o f 
theo retical s t udies for n = 3 from iKoIakowska et all (|l998l ) 
and IWinterl l|2009l ). Along with the discrete cross section 
data points, we are providing the results of our Chebyshev 
polynomial fits described in Appendix A. 
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Transfer into 3 s 




10 20 30 40 50 60 70 80 20 40 60 80 

E co i, keV E col , keV 



Figure Bl. Cross sections for excitation (left panel) and charge transfer (right panel) into 3s state. Stars show results of our calculations, 
diamonds correspond to Kolakowska (1998), and squares show results from Winter (2009). Solid line corresponds to our Chcbyshev 
polynomial fit. 




E coh keV E coh keV 

Figure B2. Cross sections for excitation (left panel) and charge transfer (right panel) into 3p state. Stars show results of our calculations, 
diamonds correspond to Kolakowska (1998), and squares show results from Winter (2009). Solid line corresponds to our Chebyshev 
polynomial fit. 




Figure B3. Cross sections for excitation (left panel) and charge transfer (right panel) into 3d state. Stars show results of our calculations, 
diamonds correspond to Kolakowska (1998), squares show results from Winter (2009). Solid line corresponds to our Chebyshev polynomial 
fit. 
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AO 


Al 


A2 


A3 


A4 


A5 


A6 


A7 


OS 


U.oll 


T fi/l 


-u.ooy 


-U. 1 to 


-u.uuoyo 


yj.loZ 


-U.UOoo 


U.UUoUo 


op 


O.OJ 


l.UO 


-U.U ( OU 






n rifi/i 

U.UU'l 


4 x 10 — ( 


-u.uuiuu 


3d 


1.59 


-0.0412 


0.0595 


-0.279 


-0.00958 


0.105 


-0.0475 


-0.0168 


4s 


-1.16 


1.08 


-0.371 


-0.0661 


-0.097 


0.191 


-0.0702 


0.015 


4p 


0.0355 


0.187 


-0.108 


-0.186 


-0.0919 


0.0756 


-0.00533 


-0.0405 


4d 


0.835 


1.4 


-0.235 


-0.0475 


-0.0222 


0.0115 


0.0116 


-0.00492 


4f 


-3.71 


-0.925 


-0.122 


-0.325 


-0.0885 


0.108 


0.016 


-0.0633 



Table Al. Fitting coefficients Aj corresponding to excitation transitions into n = 3 and n = 4 states. 





AO 


Al 


A2 


A3 


A4 


A5 


A6 


A7 


3s 


1.15 


0.963 


-1.35 


-0.242 


0.14 


0.00596 


-0.0572 


0.0564 


3p 


1.14 


-0.699 


-1.05 


-0.0922 


0.0511 


-0.0818 


-0.054 


-0.04 


3d 


-1.65 


-2.38 


-0.763 


0.137 


-0.00275 


-0.0614 


0.0162 


-0.0629 


4s 


-0.798 


1.14 


-1.24 


-0.321 


0.173 


-0.0128 


-0.0195 


-0.0199 


4p 


-2.78 


_2. 


-1.09 


-0.208 


-0.127 


-0.136 


0.0375 





4d 


-0.926 


-0.344 


-1.13 


-0.164 


0.135 


-0.0675 


-0.0372 


0.00486 


4f 


-6.95 


-3.49 


-0.617 


0.297 


0.174 


0.0395 


-0.0364 


0.000221 



Table A2. Fitting coefficients Ai corresponding to charge transfer into n = 3 and n = 4 states. 



Excitation into 4s Transfer into 4s 




10 20 30 40 50 60 70 80 10 20 30 40 50 60 70 80 

£ co i, keV £ co i, keV 



Figure B4. Cross sections for excitation (left panel) and charge transfer (right panel) into 4s state. Stars show results of our calculations. 
Solid line corresponds to our Chebyshev polynomial fit. 
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Figure B5. Cross sections for excitation (left panel) and charge transfer (right panel) into 4p state. Stars show results of our calculations. 
Solid line corresponds to our Chebyshev polynomial fit. 
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Figure B6. Cross sections for excitation (left panel) and charge transfer (right panel) into Ad state. Stars show results of our calculations. 
Solid line corresponds to our Chebyshev polynomial fit. 



Excitation into 4f Transfer into 4f 




£ co i, keV £ co i, keV 

Figure B7. Cross sections for excitation (left panel) and charge transfer (right panel) into 4/ state. Stars show results of our calculations. 
Solid line corresponds to our Chebyshev polynomial fit. 
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